Researchers have previous suggested that the response diversity of a community be measured by the diversity of responses to environmental change. For example, one can measure the response of each of the species’ intrinsic growth rate to temperature, quantify the strength and direction of these responses (e.g., as the first derivative of the response curve), and calculate the diversity of responses (e.g., by calculating variation in the first derivatives among the species in a community). When responses are nonlinear, the response diversity will be a function of the environmental state (i.e. the first derivative is a function of the value of the environmental state). So far we demonstrated this approach for quantifying response diversity in the context of a single environmental factor, but given that multiple environmental factors can change simultaneously, we need an approach that works in that context.
Owen learned about the mathematical principles from these youtube videos:
Imagine that the growth rate of a population depends on two environmental factors, e.g. temperature and salinity. We can represent the dependency as \(G = f(T, S)\), where \(G\) is growth rate, \(T\) is temperature, and \(S\) is salinity. It may be that the dependencies are linear, nonlinear, and with an interaction between temperature and salinity, hence our approach needs to be able to accommodate this phenomena.
The response of growth rate to change in temperature and salinity is the gradient / slope of this surface, with units of growth rate [per time] per temperature [degrees C] per salinity [parts per thousand]. Because the slope (first derivative) of the surface can (when dependencies are nonlinear) vary across the surface (location on the surface), and can vary in different directions on the surface, to calculate a slope we must specify the current environment (location on the surface) and the direction of change in the environment. The location on the curve is the current environmental condition, \((T_0, S_0)\), and the direction of environmental change is the unit vector \(\hat{u} = \langle U_T, U_S \rangle\).
Put another way, we calculate a directional derivative at a point on the response surface. We can write this as \(D_{\hat{u}}f(T_0, S_0)\) and can calculate it as \(f_T(T_0, S_0)U_T + f_S(T_0, S_0)U_S\), where \(f_T\) is the partial derivative of \(f(T, S)\) with respect to \(T\) and \(f_S\) is the partial derivative of \(f(T, S)\) with respect to \(S\).
Efficient evaluating in \(n\) dimensions can be done by taking the dot product of the partial derivatives at the location and the direction unit vector: \(D_{\hat{u}}f(T_0, S_0) = \triangledown f \cdot \hat{u}\) where, \(\triangledown f = \langle f_T, f_S \rangle\). (In R, the dot product of a and b is sum(a*b))
Figure 2.1 is an illustration of the principle of directional derivatives on a surface.
Figure 2.1: This figure needs considerable improvement! It is currently a keynote illustration.
Numerous mathematical functions have been used to represent how organismal performance changes with an environmental driver citation require. Moreover, multiple mathematical functions have been used to represent an interactive effect of two or more environmental drivers on species performance e.g. Thomas et al 2017. We have to decide if we are going to make simulations with different of these functions. This might increase our confidence that the method for calculating response diversity is robust to variation in types of response curve. If we decide not to, then we should likey argue that we think its robust.
Let us use the Eppley performance curve, which was used, for example, in this paper Bernhardt et al. 2018.
With one environmental variable, the performance (i.e., rate) is given by:
Adding a second environmental variable gives:
\(rate(E_1, E_2) = a_1e^{b_1E_1}(1 - (\frac{E_1 - z_1}{w_1/2})^2) + a_2e^{b_2E_2}(1 - (\frac{E_2 - z_2}{w_2/2})^2)\)
In this case, it is clear the effect of \(E_1\) and \(E_2\) is defined as being additive. For example, the value of \(E_2\) does not affect the value of \(E_1\) at which the rate is maximised (\(z_1\)), and vice-versa (see also Figure 3.1)
Figure 3.1: Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).
Including an interaction. One way to do this is to make the value of \(E_1\) at which the rate is maximised depend on the value of \(E_2\):
\(rate(E_1, E_2) = a_1e^{b_1E_1}(1 - (\frac{(E_1 + z_{int21}*E_2- z_1)}{w_1/2})^2 + a_2e^{b_2E_2}(1 - (\frac{E_2 - z_2}{w_2/2})^2\)
When \(z_{int21} = 0\) then this equation becomes the previously mentioned additive one. When \(z_{int} \neq 0\) then the value of \(E_1\) at which the rate is maximised is a function of the value of \(E_2\). We used this method for adding an interaction due to its simplicity. Other methods could be used, and if also or otherwise used could add confidence about the robustness of the method for calculating response diversity.
Figure 3.2: Nonlinear and additive dependence of a rate on two environmental variables. (a) The value of \(E_1\) at which the rate is maximised is independent of the value of \(E_2\). (b) The value of \(E_2\) at which the rate is maximised is independent of the value of \(E_1\).
First we create (or import) a table of parameter values of each species, with species in the rows and parameters in the columns. In the following example, only values of the \(z\) parameters differ among the species (which determine the location of the maximum rate).
For convenience we then convert the table of parameters into a list-column. We can then easily make performance curves of each of the species, and put those into a list-column in the same table.
## convert parameter table to a list-column of a tibble
par_list <- Partable_2_parlist(par_table)
## add performance curves
species_pars <- tibble(species = paste0("s", 1:s), pars = par_list) %>%
rowwise() %>%
mutate(expt = Make_expt(E1_series, E2_series, pars))
Here are some examples of the species’ performance curves (only with additive effects of \(E_1\) and \(E_2\)).
Figure 3.3: Performance curves for a species with maximum growth at low values of \(E_1\). Without interacting environmental effects.
Figure 3.4: Performance curves for a species with maximum growth at high values of \(E_1\). Without interacting environmental effects.
Figure 3.5: Performance curves for a species with maximum growth at low values of \(E_2\). Without interacting environmental effects.
Figure 3.6: Performance curves for a species with maximum growth at high values of \(E_2\). Without interacting environmental effects.
And now with interacting environmental effects…
For convenience we then convert the table of parameters into a list-column. We can then easily make performance curves of each of the species, and put those into a list-column in the same table.
## convert parameter table to a list-column of a tibble
par_list <- Partable_2_parlist(par_table)
## add performance curves
species_pars <- tibble(species = paste0("s", 1:s), pars = par_list) %>%
rowwise() %>%
mutate(expt = Make_expt(E1_series, E2_series, pars))
Here are some examples of the species’ performance curves (with interacting effects of \(E_1\) and \(E_2\)).
Figure 3.7: Performance curves for a species with maximum growth at low values of \(E_1\). With interacting environmental effects.
Figure 3.8: Performance curves for a species with maximum growth at high values of \(E_1\). With interacting environmental effects.
Figure 3.9: Performance curves for a species with maximum growth at low values of \(E_2\). With interacting environmental effects.
Figure 3.10: Performance curves for a species with maximum growth at high values of \(E_2\). With interacting environmental effects.
Try with and without an interaction. Therefore make two species, one with no interaction z_int = 0 and the other with z_int = 0.1. All other parameters are the same. Note that noise is added to the rate observations.
Bottom line is that the gam picks up an interaction when we have included one in the parameters used to generate the rates, and does not pick one up when we have not. This confirms that our more mechanistic thinking and methods are matching our statistical thinking and methods, and confirms that each are promising, so far.
GAM notes:
gratia package that we’ve been using to calculate derivatives of GAMs.). By the way, there is a part of this video about Model Selection, about Confidence Intervals, and about p-values for smooths (about 1h:57m to 2h:26) which I think is less useful for our purposes.method = "REML in the gam() call (the default in mgcv version 1.8-36 is method="GCV.Cp").bs = "tp")te(E1, E2) have separate marginal basis, and therefore separate smoothness parameters. They are invariant to scales of \(E1\) and \(E2\). Therefore we should be using tensor products.te(E1, E2) does not then allow inspection of the main effects and interaction term–there is essentially only one smoother which contains the main effects and interaction. So it is difficult to know if the interaction term is required. An alternative is to specify the model as s(E1) + s(E2) + ti(E1, E2), in which the third term is then only the interaction part and allows a decomposition into the different effects. Useful for examining if the interaction is required.bs = "cr".mgcvgam.check(). Only cost to large \(k\) is computational effort.Figure 3.11: Performance curves for a species without interacting environmental effects and with some noise in the rate.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rate ~ ti(E1) + ti(E2) + ti(E1, E2)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0154285 0.0004991 30.91 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## ti(E1) 3.999 4.000 14145.143 <2e-16 ***
## ti(E2) 3.937 3.997 609.233 <2e-16 ***
## ti(E1,E2) 1.001 1.002 0.023 0.883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.958 Deviance explained = 95.8%
## -REML = -5815.8 Scale est. = 0.000648 n = 2601
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rate ~ s(E1) + s(E2)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0154285 0.0003952 39.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(E1) 8.970 9.000 10199.8 <2e-16 ***
## s(E2) 7.146 8.187 475.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.974 Deviance explained = 97.4%
## -REML = -6407.3 Scale est. = 0.00040624 n = 2601
Figure 3.12: Performance curves for a species with interacting environmental effects and with some noise in the rate.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rate ~ ti(E1) + ti(E2) + ti(E1, E2)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0656082 0.0004674 140.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## ti(E1) 3.999 4.000 5449.0 <2e-16 ***
## ti(E2) 3.954 3.999 523.7 <2e-16 ***
## ti(E1,E2) 7.252 8.924 878.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.924 Deviance explained = 92.5%
## -REML = -5976 Scale est. = 0.00056817 n = 2601
##
## Family: gaussian
## Link function: identity
##
## Formula:
## rate ~ s(E1) + s(E2)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0656082 0.0009042 72.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(E1) 8.561 8.945 672.0 <2e-16 ***
## s(E2) 4.421 5.434 101.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.717 Deviance explained = 71.8%
## GCV = 0.0021379 Scale est. = 0.0021264 n = 2601
First step in calculating directional dreivatives is estimating the two partial derivatives \(f_{E1}(E1_0, E2_0)\) and \(f_{E2}(E1_0, E2_0)\) (please review the section The principle if necessary).
Partial derivative with respect to E1 (E2 is held constant)
Partial derivatives. Draw response surface for sp 1 and calculate partial derivatives at a specific location (E1 = 300, E2 = 20). To calculate the partial derivative with respect to E1, E2 must be held constant.
Figure 4.1: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.
Figure 4.2: Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20.
Partial derivative with respect to E1 when E2 is constant at 20
Figure 4.3: Partial derivative with respect to E1 when E2 is constant at 20.
Partial derivatives with respect to E2 (E1 held constant)
Figure 4.4: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.
Partial effect of E2 on the growth rate of sp 1 when E1 is held constat at E2 = 300
Figure 4.5: Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300
Partial derivative with respect to E2 when E1 is constant at 300
Figure 4.6: Partial derivative with respect to E1 when E2 is constant at 20.
Plot the two partial derivatives and relative effects
Figure 4.7: Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.
We start showing how directional derivatives can be calculated even when the direction of the environmental change is unknown. In this case, we calculate, for a specific point (E1 = 300, E2 = 20), directional derivatives in all directions.
Figure 5.1: Directional derivatives calculated in all possible direction for a specific point on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
We can measure all possible directional derivatives also for several points on the surface. This might be the case when we have multiple locations on the surface (multiple evnviromental conditions), but we do not know the direction of change.
Figure 5.2: Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
Finally, we might do the same for all the points on the surface.
Figure 5.3: Directional derivatives calculated in all possible direction for all points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
Repeating the same steps to see if it works with interactive environmental effect.
Partial derivative with respect to E1 (E2 is held constant)
Partial derivatives. Draw response surface for sp 2 and calculate partial derivatives at a specific location (E1 = 300, E2 = 20). To calculate the partial derivative with respect to E1, E2 must be held constant.
Figure 6.1: Response surface of sp1. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.
Figure 6.2: Partial effect of E1 on the growth rate of sp 2 when E2 is held constant at E2 = 20.
Partial derivative with respect to E1 when E2 is constant at 20
Figure 6.3: Partial derivative with respect to E1 when E2 is constant at 20.
Partial derivative with respect to E2 (E1 is held constant)
Figure 6.4: Response surface of sp2. The two solid lines show at which level of E1 and E2 each partial derivative is going to be calculated.
Partial effect of E2 on the growth rate of sp 2 when E1 is held constat at E2 = 300
Figure 6.5: Partial effect of E2 on the growth rate of sp 2 when E1 is held constant at E1 = 300
Partial derivative with respect to E2 when E1 is constant at 300
Figure 6.6: Partial derivative with respect to E1 when E2 is constant at 20.
Plot the two partial derivatives and relative effects
Figure 6.7: Summary plot sp1. (a) response surface of sp 1. (b) Partial effect of E1 on the growth rate of sp 1 when E2 is held constant at E2 = 20. (c) Partial derivative with respect to E1 when E2 is constant at 20. (d) Partial effect of E2 on the growth rate of sp 1 when E1 is held constant at E1 = 300. (e) Partial derivative with respect to E2 when E1 is constant at 300.
We start showing how directional derivatives can be calculated even when the direction of the environmental change is unknown. In this case, we calculate, for a specific point (E1 = 300, E2 = 20), directional derivatives in all directions.
Figure 7.1: Directional derivatives calculated in all possible direction for a specific point on the response surface of sp2. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
We can measure all possible directional derivatives also for several points on the surface. This might be the case when we have multiple locations on the surface (multiple evnviromental conditions), but we do not know the direction of change.
Figure 7.2: Directional derivatives calculated in all possible direction for several points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
Finally, we might do the same for all the points on the surface.
Figure 7.3: Directional derivatives calculated in all possible direction for all points on the response surface of sp1. Clearly, the slope of the directional derivative depends on the direction (red positive, blue negative).
Now we calculate response diversity for a community composed of 3 spp for 4 different cases: 1. Unknown direction of the environmental change 2. Direction of env change is given by the time series,and E1 and E2 change over time independently 3. Direction of env change is given by the time series,and E1 and E2 change over time with positive correlation 4. Direction of env change is given by the time series,and E1 and E2 change over time with negative correlation
[Only additives effects of the drivers on sp responses - we may change this later, so using an example with interactive effect]
Steps:
Simulate spp performance curves with the modified Eppley funciton.
Fit response surface for each sp (done with GAMs)
Data wrangling and partials derivatives calculations
Actual RD calculation
[Plotting response diversity - think about how to plot RD when we do not know the direction of env change. Surface?]
In this case the direction of the environmental change is given by the change of E1 and E2 over time.
Figure 8.1: Time series of E1 and E2 changing independently over time.
Getting partial derivatives and direction of change
Calculation response diversity
Plot response diversity over time
(ref:RD-independent-plot) Directional derivatives and response diversity with known direction of env change. E1 and E2 change independently over time. a. Species directional derivatives over time. b. Response diversity measured as similarity-based diversity metric. c. Response diversity measured as divergence (sign sensitive).
Figure 8.2: (ref:RD-independent-plot)
We first need to create a time series with E1 and E2 changing over time with negative correlation
Figure 8.3: Time series of E1 and E2 changing with negative correlation over time.
Getting partial derivatives and direction of change
Calculation response diversity
Plot response diversity over time
(ref:RD-negative-plot) Directional derivatives and response diversity with known direction of env change. E1 and E2 change with negative correlation over time. a. Species directional derivatives over time. b. Response diversity measured as similarity-based diversity metric. c. Response diversity measured as divergence (sign sensitive).
Figure 8.4: (ref:RD-negative-plot)
We first need to create a time series with E1 and E2 changing over time with positive correlation
Figure 8.5: Time series of E1 and E2 changing with positive correlation over time.
Getting partial derivatives and direction of change
Calculation response diversity
Plot response diversity over time
(ref:RD-positive-plot) Directional derivatives and response diversity with known direction of env change. E1 and E2 change with negative correlation over time. a. Species directional derivatives over time. b. Response diversity measured as similarity-based diversity metric. c. Response diversity measured as divergence (sign sensitive).
Figure 8.6: (ref:RD-positive-plot)
We use data coming from an experiment where individual ciliates species have been exposed to a gradient of nutrient, light, and their combinations in a factorial design. We first show how to calculate the partial derivatives, then we calculate the directional derivatives based on a simulated time series (in the original experiment, the level of the treatments have been kept constant throughout the expt duration). Finally, we assemble random composed communities and calculate response diversity for each of them.
(ref:ciliates-fig) Species responses to the environmental drivers. a. Species responses to nutrient concentrations. b. Species responses to light intensity.
Figure 9.1: (ref:ciliates)